function collective2

% uses 2nd ordered centered differences

clear *
clf

c=1;
gam=3;
ep=0.1;
tmaxx=2*log(3)/ep
tmax=0.992*tmaxx

%  the following parameters are adjustable
N=6000;
%M=155;

%  the following parameters are assumed fixed in the code
xR=2+sqrt(2)*c*tmax;
xL=-2-sqrt(2)*c*tmax;

% generate the points along the x-axis, x(1)=xL and x(N)=xR
x=linspace(xL,xR,N);
h=x(2)-x(1);

% generate the points along the t-axis, t(1)=0 and t(M)=tmax
k=0.99*h/(sqrt(2)*c);
lamda=c*k/h;

M=round(tmax/k)

% calculate initial values
U(1,1)=0;
V(1,1)=0;
U(N,1)=0;
V(N,1)=0;
for i=2:N-1
	U(i,1) = g0(x(i));
	V(i,1) = h0(x(i));
	f0 = gam*g0(x(i))*h0(x(i)) - h0(x(i));
	dh = (h0(x(i+1)) - h0(x(i-1)))/(2*h);
	dg = (g0(x(i+1)) - g0(x(i-1)))/(2*h);
	ddh = (h0(x(i+1)) - 2*h0(x(i)) + h0(x(i-1)))/h^2;
	ddg = (g0(x(i+1)) - 2*g0(x(i)) + g0(x(i-1)))/h^2;
	U(i,2) = g0(x(i)) - 2*c*k*dh + k^2*(c^2*ddg - c*ep*(gam*dg*h0(x(i)) + gam*g0(x(i))*dh - dh));
	ccv = -2*gam*c*h0(x(i))*dh + gam*g0(x(i))*(-c*dg + ep*f0) + c*dg - ep*f0;
	V(i,2) = h0(x(i)) + k*(-c*dg + ep*f0) + k^2*(c^2*ddh + 0.5*ep*ccv); 
end;

%plot(x,U(:,2),'-r')
%pause

for j=3:M
	U(1,j)=0; V(1,j)=0;
	U(N,j)=0; V(N,j)=0;
	for i=2:N-1
		vv1 = -2*gam*lamda*V(i,j-1)*(V(i+1,j-1) - V(i-1,j-1)) - V(i,j-2)*(gam*U(i,j-1) - 1);
		vv2 = 2*lamda^2*V(i+1,j-1) + 2*(1-2*lamda^2)*V(i,j-1) + 2*lamda^2*V(i-1,j-1) - V(i,j-2);
		vv3 = 1 - 0.5*ep*k*(gam*U(i,j-1) - 1);
	  	V(i,j) = (vv2 + 0.5*ep*k*vv1)/vv3;
		U(i,j) = U(i,j-2) - 2*lamda*(V(i+1,j-1) - V(i-1,j-1));
	end;
end;

%  calculate asy solution
t=(M-1)*k;
for i=1:N
	theta1 = x(i) - sqrt(2)*c*t;
	theta2 = x(i) + sqrt(2)*c*t;
	di1 = g0(theta1) + sqrt(2)*h0(theta1);
	di2 = g0(theta2) - sqrt(2)*h0(theta2);
	et = exp(ep*t/(1+c));
	gf = gam*(1-et);
	ue(i) = di1/(di1*gf+2*et) + di2/(di2*gf+2*et);
	ve(i) = (di1/(di1*gf+2*et) - di2/(di2*gf+2*et))/sqrt(2);
end;

% plot results
hold on
box on
plot(x,U(:,M),'-r')
plot(x,V(:,M),'b')
plot(x,ue,'--r')
plot(x,ve,'--b')
%axis([xA xB -1 1])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);
%legend([' Upwind: M = ', num2str(M)],' Exact',4);
set(findobj(gcf,'tag','legend'),'FontSize',12,'FontWeight','bold'); 
hold off



% initial g(x)
function q=g0(x)
% step between -0.5 to 0.5
%q=0.5*(sign(x+0.5) + 1) - 0.5*(sign(x-0.5) + 1);
q=0;
if abs(x)<1 
	q=0.5*(cos(pi*x)+1);
end


% initial h(x)
function q=h0(x)
q=0;
















